c----------------------------------------------------------------------
      real*8 function axb_vec (a, b, n)
      integer i, n
      real*8 a(n), b(n)
      axb_vec = 0d0
      do i = 1, n
         axb_vec = axb_vec + a(i) * b(i)
      end do
      return
      end
c-----------------------------------------------------------------
      subroutine dc_axb (a, b, c, itr, na, nb, nc, na1, na2, nb2)
      real*8 a(na, na2), b(nb, nb2), c(nc, nb2)
      integer na, nb, nc, na1, na2, nb2, i1, i2, j, itr
c     C = A'* B  if itr.eq.1
c     C = A * B  if itr.ne.1
      if (itr .eq. 1) then
         do i1 = 1, na2
            do i2 = 1, nb2
               c(i1, i2) = 0d0
               do j = 1, na1
                  c(i1, i2) = c(i1, i2) + a(j, i1) * b(j, i2)
               end do
            end do 
         end do
      else
         do i1 = 1, na1
            do i2 = 1, nb2
               c(i1, i2) = 0d0
               do j = 1, na2
                  c(i1, i2) = c(i1, i2) + a(i1, j) * b(j, i2)
               end do
            end do 
         end do
      end if 
      return
      end
c-----------------------------------------------------------------
      subroutine da_bxv (a, b, v, itr, na, nb, na1, na2, nb2)
      real*8 a(na, na2), b(nb, nb2), v
      integer na, nb, na1, na2, nb2, i1, i2, j, itr
c     A = B'* v  if itr.eq.1
c     A = B * v  if itr.ne.1
      if(itr .eq. 1) then
         if(nb2 .ne. na1) then
            print*, 'da_bxv-error, please check na1, nb2', na1, nb2
            stop
         end if 
         do i1 = 1, na1
            do i2 = 1, na2
               a(i1, i2) = b(i2, i1) * v
            end do 
         end do
      else
         if(nb2 .ne. na2) then
            print*, 'da_bxv-error, please check na2, nb2', na2, nb2
            stop
         end if 
         do i1 = 1, na1
            do i2 = 1, na2
               a(i1, i2) = b(i1, i2) * v
            end do 
         end do
      end if 
      return
      end
